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ABSTRACT 

HI intensity mapping is an emerging tool to probe dark energy. Observations of the redshifted 
HI signal will be contaminated by instrumental noise, atmospheric and Galactic foregrounds. 
The latter is expected to be four orders of magnitude brighter than the HI emission we wish 
to detect. We present a simulation of single-dish observations including an instrumental noise 
model with 1// and white noise, and sky emission with a diffuse Galactic foreground and 
HI emission. We consider two foreground cleaning methods: spectral parametric fitting and 
principal component analysis. For a smooth frequency spectrum of the foreground and in¬ 
strumental effects, we find that the parametric fitting method provides residuals that are still 
contaminated by foreground and 1// noise, but the principal component analysis can remove 
this contamination down to the thermal noise level. This method is robust for a range of dif¬ 
ferent models of foreground and noise, and so constitutes a promising way to recover the 
HI signal from the data. However, it induces a leakage of the cosmological signal into the 
subtracted foreground of around 5%. The efficiency of the component separation methods 
depends heavily on the smoothness of the frequency spectrum of the foreground and the 1// 
noise. We find that as, long as the spectral variations over the band are slow compared to the 
channel width, the foreground cleaning method still works. 

Key words: diffuse radiation - large-scale structure of Universe - cosmology: observations 
- methods: statistical - radio continuum: general, galaxies - radio lines: galaxies, ISM 


1 INTRODUCTION 

One of the main challenges of modern cosmology is to explain 
the late-time acceleration of the expansion of the Universe. There 
are several independent methods to probe dark energy: Baryonic 
Acoustic Oscillations (BAOs), weak and strong gravitational lens- 
ing, cluster counts and supernovae. However, BAO measurements 
appear to be the most powerful cosmological tool at low redshift be¬ 
cause they are limited by statistical rather than systematics errors. 
Using the BAOs as a standard ruler measures the expansion of the 
Universe as a function of redshift, and so constrains the properties 
of dark energy (e.g. Weinberg et al. 2013). 

A complementary method to the usual large optical surveys of 
galaxies to study BAOs is HI intensity mapping (Battye, Davies & 
Weller 2004; Peterson, Bandura & Pen 2006; Chang et al. 2008; 
Loeb & Wyithe 2008; Peterson et al. 2009). This method aims to 
give a tomographic distribution of the neutral HI emission present 
in the recent Universe over large angular scales. Simulations indi¬ 
cate that the HI intensity mapping technique will give very precise 
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constraints on cosmological parameters, and in particular, on the 
dark energy equation of state at low redshift (Chang et al. 2008), 
and at high redshift (McQuinn et al. 2006; Bowman, Morales & 
Hewitt 2007; Mao et al. 2008). This sensitivity comes from the 
large volume of the survey. 

Some HI experiments are currently underway such as BAO- 
radio', BINGO^, CHIME^, FAST^, TIANLAr. Using intensity 
mapping, the Green Bank Telescope (GBT)® has provided the first 
detection of the HI signal at z ~ 0.8 cross-correlated with the Wig- 
gleZ Dark Energy Survey (Masui et al. 2013). This detection shows 
the HI intensity mapping as a promising tool and gives a lower limit 
on the fluctuations power of the HI signal. The first phase of the 
SKA instrument will be built in the next decade and it will offer 
a broad range of frequencies and a large survey area. Thus, this 
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instrument has great potential to deliver maps of the HI intensity 
(Bull et al. 2015). 

To convert the promise into reality radio observations will 
have to deal with different contaminants, which can dominate the 
signal of interest in the data, such as astrophysical foregrounds, ra¬ 
dio frequency interference (RFI) and instrumental noise. Another 
important contaminant is time variable noise introduced during 
propagation of the signal through the atmosphere, which gives an 
additional contribution to the 1// noise of the instrument. The am¬ 
plitude of the atmospheric effects depends on the observing fre¬ 
quency, on the elevation of the instrument above sea level and on 
the instantaneous weather conditions. 

The most important challenge for any intensity mapping ex¬ 
periment is the control of foreground emissions. Thus, the data 
analysis must include a robust foreground subtraction algorithm. 
At ~ 1 GHz, the most relevant foregrounds are a combination of 
Galactic emission, mostly synchrotron, and that from the back¬ 
ground of extragalactic point sources. See, for example, the discus¬ 
sion in Battye et al. (2013). These emissions are at least four orders 
of magnitude larger (7^ 1000 mK) than the HI signal fluctuations 

(dTb 1 mK). In order to subtract the foreground, the high spectral 
resolution offered by any HI experiment allows us to exploit the fre¬ 
quency information. In particular the spectra of the foregrounds are 
expected to be smooth and can be approximated over the frequency 
range of interest to first order by a modified power-law with a spa¬ 
tial curvature as a function of frequency (Kogut 2012). This spectral 
smoothness can be used to separate the HI signal from foreground 
signals. The most common approach is to fit a smooth function to 
the data in the frequency space and remove it. 

Several component separation techniques have been discussed 
in the literature for removal of the Galactic foreground, both spe¬ 
cific to low redshift intensity mapping and to epoch of reionisa¬ 
tion experiments. We can classify the foreground cleaning methods 
into two different categories: parametric and blind methods. The 
parametric methods are model-dependent and consist of applying 
a parametric fitting to each pixel in the maps (Ansari et al. 2012). 
The blind methods do not require any assumption on the properties 
of the foregrounds and the instrument response. Examples of such 
methods are FASTICA (Chapman et al. 2012; Wolz et al. 2014), 
Correlated Component Analysis (CCA) method (Bonaldi & Brown 
2014), Karhunen-Loeve Decomposition (Shaw et al. 2014), GMCA 
(Chapman et al. 2012), principal component analysis (PCA) and 
independent component analysis (ICA) (Alonso et al. 2015). All 
methods must also deal with instrumentally induce effects such as 
mode-mixing of angular and frequency fluctuations induced by the 
dependence of the beam with frequency. The component separation 
methods are based on the spectral smoothness of the foregrounds, 
so the calibration will be a critical step in order not to compromise 
this smoothness. 

There are two instrumental approaches to HI mapping, ei¬ 
ther using single dishes with multiple feeds or interferometer ar¬ 
rays. The single-dish approach offers a relatively cheap and sim¬ 
ple way for doing intensity mapping. Unlike single-dish exper¬ 
iments which naturally have good surface brightness sensitivity, 
interferometer arrays suitable for intensity mapping will require 
many close-packed elements in order to detect the very low sur¬ 
face brightness HI signal and thus need big correlators (Thomp¬ 
son, Moran & Swenson 2008). Both kinds of experiment will have 
to deal with potential systematics similar to those encountered in 
cosmic microwave background (CMB) imaging experiments, such 
as gain variations, spillover, ground pickup and correlated noise in 
frequency. Single-dish experiments will require stable receiver sys¬ 


tems, good calibration and an appropriate scanning strategy. On the 
other hand, interferometers are known to deal more naturally with 
systematics and with foregrounds than single dishes and hence re¬ 
ceiver stability (e.g. Dickinson et al. 2004; Readhead et al. 2004), 
etc. is not such an important issue. However, we point out that 
existing interferometers are limited by the small number of their 
smallest baselines and hence fail to provide the required surface 
brightness sensitivity (Bull et al. 2015). 

In this paper, we focus on the concept of a single-dish ex¬ 
periment with the BINGO (BAO from Integrated Neutral gas Ob¬ 
servations) instrument, which aims at mapping the HI emission at 
frequencies from 960 to 1260 MHz (z = 0.12 — 0.48). This exper¬ 
iment will measure the HI power spectrum, and will detect for the 
first time the B ary on Acoustic Oscillations around 1 GHz. Some of 
the details of the BINGO experiment can be found in Battye et al. 
(2013). Though we use the BINGO instrument as a concrete ex¬ 
ample, nearly all the following analysis concepts can be applied to 
other single-dish instruments. 

We organise the paper as follows. In Section 2, we describe 
our simulations in which we incorporate foreground and instru¬ 
mental noise models. We also make some predictions of the to¬ 
tal atmospheric contribution to the instrumental noise level and of 
the atmospheric fluctuations coming from the inhomogeneous dis¬ 
tribution of water vapour. In Section 3, we focus on two simple 
foreground and noise subtraction procedures: parametric fitting and 
principal component analysis. The success of these methods de¬ 
pends on the smoothness of the frequency spectrum of the noise 
and the foreground. In Section 4, we place some requirements on 
the smoothness of these spectra needed for extracting the HI signal 
from the instrumental l/f noise and from the brighter foreground, 
the synchrotron emission. In this way, we assess the robustness of 
the cleaning methods, according to the smoothness of the frequency 
spectrum of the simulated data. 


2 SIMULATION OF A SINGLE-DISH EXPERIMENT 

In order to explore the foreground cleaning methods we use simu¬ 
lations of the proposed BINGO telescope (Battye et al. 2013) as a 
concrete example of a single-dish instrument. To do this we require 
simulated maps of the sky at the observation frequencies. We pro¬ 
duce a time-ordered data stream with a foreground model which in¬ 
cludes Galactic synchrotron plus a background of unresolved point 
sources, detailed in Section 2.1, while the way we produce the HI 
signal is described in Section 2.2. In Section 2.3, we introduce the 
map-making method used to obtain the sky maps of the experiment 
and describe the instrumental noise model (1// and thermal noise) 
in Section 2.4. Finally, in Section 2.5, we make some predictions 
of the amplitude of the atmospheric noise. 

2.1 Foreground model 

2.1.1 Galactic synchrotron 

To generate a template of the sky emission, we use the reprocessed 
408 MHz Haslam et al. map (Remazeilles et al. 2014), which con¬ 
stitutes a good tracer of the diffuse synchrotron emission. The syn¬ 
chrotron spectrum in terms of brightness temperature can be ap¬ 
proximated by T(y) oc vP+<^ln(v/vo) (Kogut 2012), where v is the 
radiation frequency, C the curvature deflned with respect to a ref¬ 
erence frequency Vq and P is the spectral index at v = Vq. Obser¬ 
vations have indicated that there are spatial variations of the spec- 
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Table 1. Summary of the different models of the Galactic synchrotron emis¬ 
sion. 


Synchrotron 

Characteristics 

Mean p 

r.m.s. p 

Model 1 

P constant on the sky 

-2.8 

P constant 

Model 2 

A Gaussian spatial distribution of p 

-2.8 

0.17 

Model 3 

de Oliveira-Costa et al. (2008) model 

-2.5 

0.03 


tral index (Reich & Reich 1988; Platania et al. 1998; Davies et al. 
2006). We extrapolate this template to the frequencies of interest 
by using 3 different models for P listed below, from the simplest to 
the more complicated ones: 

1. We ignore any variation across the sky of the spectral index p, 
we fix this index to an average value estimated at frequencies near 
- 1 GHz p = -2.8 (Platania et al. 1998). 

2. We assume a Gaussian spatial distribution of the synchrotron in¬ 
dex p, with <p>=-2.8 and a r.m.s. value of 0.17 (Platania et al. 
1998). 

3. We use spectral model of the global sky from 10 MHz to 
100 GHz developed by de Oliveira-Costa et al. (2008). This final 
model is the most realistic one. It includes spatial correlations of the 
Galactic emission across the sky and a frequency curvature modifi¬ 
cation of the synchrotron index p. The mean value of P is —2.5 and 
the steepening of this index is 0.03. 

The models listed above are summarised in Table 1 . 

A high-resolution template of synchrotron emission is re¬ 
quired to make realistic tests of foreground removal methods. How¬ 
ever, the resulting synchrotron map has a resolution correspond¬ 
ing to a beam with FHWM equal to 56arcmin. We require a 
higher resolution map so it cannot be directly used as a template 
of the synchrotron emission. And, there are no other full-sky as¬ 
tronomical data sets with resolution better than ~ 1°. Hence, to 
account for these small-scale fluctuations, we add to the origi¬ 
nal map a random Gaussian realisation with a power spectrum 
Q = £7(exp(—where y = —2.7, Gsim is the Gaussian 
width of the simulation and £ the multipole. The details are given 
in Miville-Deschenes et al. 2007 and Remazeilles et al. 2014. 


2.7.2 Extragalactic point sources 

We assume that the distribution of such sources is not spatially 
correlated, that is to say the clustering is weak (Liu, Tegmark & 
Zaldarriaga 2009) and hence that they are Poisson distributed. The 
clustering part increases the pixel-pixel correlations (Battye et al. 
2013) and thus can have an impact on the foreground removal 
method. In subsequent work, we will investigate this contribution. 
Extragalactic point sources can be divided into two populations. 
The first component comprises bright and isolated point sources 
that can be readily detected by the instrument and removed directly 
using the data of the experiment. The second population consists of 
a continuum of unresolved sources. 

At radio frequencies (GHz), the r.m.s. confusion Gc in a tele¬ 
scope beam with the full width at half maximum 0fwhm can be 
approximated by (Condon 1974) 



For the BINGO telescope, with a 0fwhm = 40arcmin, this is 
around 320 mJy at 1000 MHz, thus BINGO will be subject to con¬ 


fusion noise when considering a continuum detection. This is irrel¬ 
evant for an HI line signal. 

The brightness of each source is drawn from the differential 
source counts ^, with the number of sources per steradian N and 
per unit of flux S. In Battye et al. (2013), they use data from multi¬ 
ple continuum surveys at 1.4 GHz (Mitchell & Condon 1985; White 
et al. 1997; Ciliegi et al. 1999; Gruppioni et al. 1999; Hopkins 
et al. 1999; Richards 2000; Bondi et al. 2003; Fomalont et al. 2006; 
Owen & Morrison 2008; Seymour et al. 2008; Ibar et al. 2010) and 
fit a 5th order polynomial to these data 


logic 


) 




( 2 ) 


where oq = 2.593, ai = 9.333 x 10 02 = —4.839 x 10 = 

2.488 X 10-^ a 4 = 8.995 x 10-^ andas = 8.506 x 10-^; and Aq = 
1 Jy^/^ sr“^ and Sq = I Jy. The power-law spectral function with a 
Gaussian distributed index is given by 


S(v).S(l,40H.)(^)-“ (3) 

The spectral index a is randomly chosen from a Gaussian distribu¬ 
tion 


P{a) 




=exp - 


(a-ao)^l 
2o2 J ’ 


(4) 


with a mean of Gq = 2.5 and a width distribution of G = 0.5 
(Tegmark et al. 2000). 

Assuming that the sources with flux S > iSmax can be sub¬ 
tracted from the data, we estimate the mean brightness temperature, 
contributed due to the remaining sources, by 


(v, w) - j^ ^pix 11^ (y) ’ 


(5) 


where St is the flux of the point source i at 1.4 GHz and Ilsky the 
pixel size qual to 0.22 arcmin^. The parameter dB/dT = is 

the conversion factor between intensity units to brightness temper¬ 
ature units, being the Boltzmann constant, X the wavelength of 
the incoming radiation and ^Ipix the pixel size. 

We expect that the brightest sources could be removed directly 
from the BINGO data or could be masked using the NR AO VLA 
Sky-Survey (NVSS) (Condon et al. 1998), which is considered to 
be 99% complete at a flux density limit of 3.4 mJy. For our sim¬ 
ulation, to be conservative, we take iSmax = lOOmJy, which corre¬ 
sponds to ~ 1 source per square degree. We expect to either subtract 
or mask most of the brightest radio sources above this flux density. 
We will investigate in a following paper the residual contribution 
due to the variability of radio sources and calibration issues. 

In the following, the maps are created using the HEALPix pix- 
elisation scheme (Gorski et al. 2005). A foreground map of this 
simulation at 1000 MHz is given in Fig. 1. The colour bar repre¬ 
sents the brightness temperature in mK. 


2.2 HI signal 

We use the CORA software (Shaw et al. 2014) to simulate the HI 
brightness temperature. We assume the Planck-i-WP-i-highL-i-BAO 
cosmological model given in Planck Collaboration et al. (2014). 
The HI brightness temperature can be written as a sum of two parts 

rb = fb(i+5Hi). (6) 
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Table 2. Instrumental parameters for BINGO simulation. 


Survey parameters 

Redshift range [Zmin,Zmax] 

[0.13,0.48] 

Frequency range [VminAmax] (MHz) 

[960, 1260] 

Channel width Av (MHz) 

15 

FWHM (arcmin) at 1 GHz 

40 

Number of feed horns nf 

56 

Sky coverage flsur (deg^) 

3000 

Observation time fobs (yr) 

1 

System temperature T^ys (K) 

50 

Sampling rate (Hz) 

0.1 


Figure 1. Mollweide projection of the foreground with a background 
of unresolved point sources (S < lOOmJy) and synchrotron emission at 
1000 MHz in celestial (RA/Dec) coordinates with RA=0° at the centre and 
increasing to the left. The white solid lines define the spaces expected to be 
observed by the BINGO experiment. 


where 5 hi is the HI density contrast and Tj, mean HI brightness 
temperature given by 


fb(z)=0.3K 


^HI \ / ^m + (l+z) ^ 

TF3 J V 0^9 J / 


(7) 


We assume that the neutral HI fraction is ^2 hi = 5 x 10 ^ (Switzer 
et al. 2013) and the HI bias is independent of scale and redshift with 

bill = 1- 

The HI brightness temperature power spectrum can be mod¬ 
eled as 


Pt, {k,z) = {z)[bm + fp^fD^iz)P„,ik,z), (8) 


where k\\/k with the flat-sky approximation, P^nik^z) the mat¬ 
ter power spectrum, D(z) the linear growth factor normalised by 
T)(0) = 1, and / the linear growth rate / = dlogD/dloga, where 
a is the cosmological scale factor. The HI angular power spectrum 
is obtained from Gaussian random fields with the flat sky angular 
power spectrum (Datta, Choudhury & Bharadwaj 2007) 


cf'(Av) = rd^||Cos(A:||rvAv)Pr,(k). (9) 

Ttr Y Jo 


where is the comoving distance, k has components k\\ and i/ry 
along the line-of-sight and in the plane of the sky respectively. Us¬ 
ing these inputs, we generate the maps of the HI signal which have 
r.m.s. fluctuations around 0.1 mK. 


2.3 Simulation of a single-dish experiment 

We consider a single-dish experiment based on the BINGO con¬ 
cept. BINGO will be a dual mirror Compact Antenna Test Range 
(CATR) telescope with a 40 m primary mirror and an offset focus. 
Apart from the telescope optics the design of the instrument is simi¬ 
lar to that of Battye et al. (2013). The proposed BINGO experiment 
will have a receiver array containing between 50 and 60 feed horns. 
In our simulation, we model the receiver plane with 56 feed horns 
with a 90 m focal length. We consider the frequency range from 
960 MHz (z = 0.48) to 1260 MHz (z = 0.13). To decrease the com¬ 
putational speed, we choose to divide the 300 MHz band into 20 
channels, each of 15 MHz bandwidth, though the actual instrument 
will have much narrower frequency channels to facilitate RFI ex¬ 
cision. The sampling rate is 0.1 Hz. The instrumental parameters 
used for our simulation are listed in Table 2. 


We assume that the horns are arranged in a rectangular con¬ 
figuration spaced 3.3 m apart and the beams are given by a circular 
Gaussian. The beams are diffraction-limited and, therefore, the full 
width at half maximum OpwHM of the beam can be scaled to any 
frequency V by 

Vn 

0FWHm(v) = 0FWHm(Vo) —, (10) 

with Vo = 1000 MHz and 0fwhm(vo) = 40 arcmin. 

For the following simulations, we will assume that the tele¬ 
scope will map a 15° declination strip centred at —5° as the sky 
drifts past the telescope. The declination of —5° has been chosen 
to minimise the foreground emission, which is lowest between 10 
and —10° declination. We assume one full year of on-source inte¬ 
gration. In practice, this will likely represent about 2 years of real 
observation time since we could consider only night observations 
we will probably remove some data due to foreseeable technical 
issues like such as radio frequency interference, weather downtime 
etc. 

In order to obtain the simulated BINGO maps, we use a 
maximum likelihood map-making algorithm (Stompor et al. 2002; 
Hamilton 2003). We model the timelines d as d = As + n, where s 
is the pixelized sky signal which is mapped into the timelines and 
corrupted by noise n. The pointing information is represented by 
the pointing matrix A of size Nsamples x Npixeis. which connects 
the time index to the pixel index. The map-making step is given by 

s = {A^N~^A)~^A^N~^d, (11) 

where N is the noise covariance matrix and s is the best estimate 
of s. An impact of the l/f noise is to induce slow drifts of the 
gains of the receivers. If we do not take steps to mitigate it, the 1// 
noise will introduce stripes in the maps along the direction of the 
drift scan. The inversion of (A^A^“^A) is performed by using the 
preconditioned conjugate gradient method. The preconditioner is a 
pixel domain diagonal matrix weighting the pixels by the number of 
times they have been observed. This method is described in detail 
in Cantalupo et al. (2010). 

We set the HEALPix resolution of the map equal to nside = 
128, which corresponds to a map pixel size of 27 arcmin. The fo¬ 
cal plane configuration will lead to some gaps in the observed sky 
band. To correct for this, we rotate the beams of the horns on the 
sky with an angle ~ 5°. In Fig. 2, we show the drift scan strips of 
the sky emission. In the following, we consider a single frequency 
channel centered at 997.5 MHz to display the results. The top panel 
shows the HI signal and the bottom panel the Galactic synchrotron 
emission plus a background of unresolved point sources. The am¬ 
plitude of the foreground emission is much higher than the sig¬ 
nal of interest by four orders of magnitude (note the difference in 
colour-scale in the strips). These maps are plotted with a Cartesian 
projection using the HEALPix software. 
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Figure 2. Two drift scan strips of the observed sky with only HI signal in 
the top panel while the bottom panel includes Galactic synchrotron emis¬ 
sion plus a background of unresolved point sources. The mask of the Galac¬ 
tic plane at \b\ < 20° is plotted with grey pixels. Note the different linear 
intensity scales. Colour bars represent the temperature in mK. 


2.4 Instrumental noise 


The timelines are corrupted by thermal (white) noise. The optimal 
sensitivity of the BINGO experiment per pixel can be defined as 
follows 


at = 


^sys 

\/^pix^^ 


( 12 ) 


where Av is the frequency channel given in Table 2. We assume the 
same system temperature for all receivers. The parameter 
is the integration time per pixel defined by 

^pix = ^f^obs ^ 5 (13) 

*^sur 


where Uf denotes the number of feed horns, ?obs is the total inte¬ 
gration time, Osur is the survey area and Opix is the beam area. 
The values of these parameters are given in Table 2. We assume 
^pix = integration time of one year, we obtain 

at = 25 a/K. 

Our simulation also contains 1// noise, produced by gain fiuc- 
tuations of the amplifiers and thus this noise is correlated across 
all frequency channels. The impact of these fiuctuations is usually 
simulated in the frequency domain using a 1// power spectrum. To 
generate realistic 1// noise we use an algorithm to create a time 
sequence of white noise. Then we compute its Fourier Transform 
and weight the data by a zero-mean power spectral density with the 
distribution 


Ad — : 


samp 


1 + 




(14) 


where Vsamp is the sampling frequency. / is the discrete Fourier 
transform sample frequency, given a length of the number of time 
samples and a sample spacing of 1/Vsamp- The 1// knee value is 
the integration time (the inverse of the discrete Fourier transform 
sample frequency) at which the thermal and 1// noise make equal 
contributions to the power spectral density. Finally, we compute the 
inverse Fourier transform of these data to obtain the time-ordered 
data of the noise. 

In practice, we will filter the data on timescales of 20 min 
during data processing; this is the timescale for which the largest 
structures of interest take to drift through the BINGO field-of-view. 
For simplicity, we simulate each 20 min timestream separately and 
join them together by fixing the first sample of the nth equal to the 
(n — 1) sample. We assume the same value of the knee frequency 
for each receiver /^nee = Hz for a 15 MHz channel bandwidth, 
which corresponds to the value we aim to achieve with the BINGO 
pseudo-correlation receivers. The 1// slope index a is assumed to 
be 1. In this paper, we start by making the assumption that the 1// 
noise is perfectly correlated between the frequency channels, which 
is what is expected if we are dealing with simple gain fiuctuations 


and we assume a fiat frequency spectrum (P = 0). In Section 4, we 
investigate what happens if these assumptions are relaxed. 

In Fig. 3, we show two maps of the noise with thermal noise 
only (in the top panel) and with added 1// noise (in the bottom 
panel). One can notice the stripes along the direction of the scan 
induced by the instrumental 1// noise, which is much larger than 
the thermal noise by a factor of ~ 100. 


2.5 Atmospheric noise 

The observations of a ground-based telescope are affected by the at¬ 
mosphere at different levels, depending on the observing frequency 
and of course on the prevailing weather conditions. The incoming 
signal will be absorbed or scattered by the atmosphere. This effect 
increases the noise level of the instrument and we quantify this in 
Appendix Al. Around 1 GHz, the optical depth is dominated by 
oxygen which is constant in time with a small contribution from 
water, usually quantified in terms of the precipitable water vapour 
(PWV). Hence, it can undergo changes on an hourly (or faster) 
timescales and vary spatially. In Appendix Al we show that the 
amount of PWV at 1 GHz does not have a significant impact on the 
total brightness temperature of the atmosphere. We find an atmo¬ 
spheric contribution ~ 1.81 K in days with favourable atmospheric 
conditions PWV < 2 mm, and ^ 1.82 K in days with unfavourable 
weather PWV > 4 mm. Thus, we expect the atmospheric contribu¬ 
tion to the system temperature to be at the level of a few percent. 

However, we are also concerned with the fluctuating part of 
the emission and absorption, which leads to an additional 1// noise¬ 
like component in the time-ordered data. We quantify this effect in 
Appendix A2. The main source of fluctuations is atmospheric tur¬ 
bulence in the troposphere. Oxygen molecules are uniformly mixed 
in the atmosphere, but the water vapour molecules show inhomo¬ 
geneities. Even if the water contributes a small percentage to the 
noise level, around 1%, its variations with time can produce fluctu¬ 
ations that are large compared to the HI signal and thermal noise. 
This can generate some correlations between the time streams of 
data from different receivers and for the same receiver at differ¬ 
ent times. These fluctuations depend on the weather conditions at 
the observing site, on the telescope, on the frequency bands and 
on the observing strategy (Church 1995; Lay & Halverson 2000). 
The effects of atmospheric noise have been well studied in cos¬ 
mic microwave background experiments (e.g. Davies et al. 1996; 
Sayers et al. 2010), and can be approximated by 1// noise at low- 
frequency. However, the amplitude of the atmospheric fluctuations 
in the time streams has not been measured at 1 GHz. From the 
back-of-the-envelope calculation detailed in Appendix A2, we find 
that the amplitude of the atmospheric fluctuations in antenna tem¬ 
perature for a single dish is ATatm ^0.01 mK. Thus, the ampli¬ 
tude of these fluctuations is below the instrumental noise, which 
is expected to be ~1 mK for a frequency resolution of 15 MHz for 
BINGO. So, the contribution from the atmospheric fluctuations ap¬ 
pears not to be a challenge for a single-dish experiment observing 
around 1 GHz when the weather is stable. 


3 FOREGROUND AND INSTRUMENTAL NOISE 
SUBTRACTION 

For an observing frequency around 1 GHz, the synchrotron emis¬ 
sion and the extragalactic point sources are the most relevant fore¬ 
grounds. The removal of the foregrounds and instrumental 1// will 
rely on the smoothness of their frequency spectra. In this section. 
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Figure 3. Maps of the drift scan of the instrumental noise. The top panel represents the thermal noise and the bottom panel the 1 // noise. Note the different 
intensity scales. Colour bars represent the temperature in mK. 


we want to quantify how well the foregrounds can be subtracted 
in the presence of thermal and 1// noise. Our philosophy is to fo¬ 
cus on two simple cleaning procedures, parametric fitting and a 
blind method with principal component analysis (PCA). We de¬ 
scribe these methods in Section 3.1 and present their results in 
Section 3.3. We also demonstrate the possibility of using a blind 
method to remove the instrumental 1// noise in Section 3.2. In the 
following, we assume no systematics and a perfect calibration of 
the data. 


Cij = = wi T(yi^p)T{yj,flp). (17) 

p=\ 

where Tiyi^hp) is the brightness temperature along the direction of 
the line-of-sight hp and for the frequency channel V/. Therefore, we 
can compute the entries of the correlation matrix between each pair 
of frequency channels 


3.1 Methods 


3.1.1 Parametric fitting 


Parametric fitting is a common method to parameterise foregrounds 
(e.g. Brandt et al. 1994; Ansari et al. 2012). The approach of the 
method is to fit directly an explicit parametric model of the fore¬ 
grounds and noise to each pixel of the maps along the frequency 
direction. The common foreground model is a modified power-law. 
As the main foreground emission, the Galactic synchrotron can be 
approximated by a parametric distribution with a curvature to first 
order (Kogut 2012). 

The i-th pixel of the simulated map of the sky at the frequency 
Vj can be written as the sum of the intensity of the HI signal ^21 cm’ 
the foreground emissions and the noise of the instrument 

fh J + %J + Kj- ( 15 ) 

The hat symbol denotes a modelled quantity. We make the assump¬ 
tion that the foreground ffg and l/f noise can be modelled by 


f4^. + f^-+=AM - 


nj 


Vo 


(16) 


where P is the spectral index and A is the amplitude in mK. This 
assumption on the spectral slope of the l/f noise can be justified 
by the fact that the l/f noise fluctuations are expected to have a 
spectral form similar to the system temperature, which can be ap¬ 
proximated by a power-law over the BINGO frequency range. We 
fit Eq. 16 for each pixel of the map in the frequency direction min¬ 
imised using a least-squares method. 


3.1.2 Principal Component Analysis (PC A) 

PC A (Murtagh & Heck 1987) has the advantage of being a non- 
parametric method and so requires no specific prior information on 
the spectra of the foreground and the noise. This method consists of 
transforming the independent maps of each frequency channel into 
orthogonal modes according to the covariance between frequen¬ 
cies. We consider the data to be a matrix 5, with Nf x Np elements. 
Nf denotes the number of frequency channels and Np the number 
of pixels in the map. We compute the frequency covariance matrix 
from the simulated data 


^jj ^kk 

where the indices run from 1 to Nf. We diagonalise the correlation 
matrix of the full data set with an eigenvalue decomposition and 
obtain 


P"i?P = A = diag{)Li,....)LivJ, (19) 

where the diagonal elements of the matrix A are the eigenvalues 
Xj of the matrix R and the matrix P is an orthogonal matrix which 
contains the eigenvectors. The variance of each mode is given by 
the amplitude of the eigenvalues Xj, so each eigenvalue measures 
the contribution of its corresponding eigenvector to the total sky 
variance. 

This method parameterises the foreground and noise compo¬ 
nents and produces independent eigenfunctions, which convert the 
spectral correlation into a number of largest variance modes. We 
pick the eigenvalues with the correlated components in frequency 
with the larger variances. So, we build a matrix P^ with only the 
corresponding eigenvectors and we use this matrix to decompose 
the data into eigenfunctions (|) 

(^ = P/S. (20) 

The maps Sc of the reconstructed foreground and l/f noise 
are obtained by transforming back to the frequency space 

Sc = Pc^. (21) 

Finally, we find the maps of the reconstructed HI signal S'hi by 
subtracting the input maps and the reconstructed foreground and 
l/f noise 

Sui = S-Sc. ( 22 ) 


3.2 l/f noise subtraction using PC A 

First, we apply the PC A method to thermal and l/f noise com¬ 
ponents only, ignoring foregrounds for the moment. The frequency 
spectrum of the correlated noise (1//) is also expected to be smooth 
in frequency. Thus, one can use the PCA method to remove the 
instrumental l/f noise. We show the result in Section 3.2.1 and 
we test the robustness of this noise removal method with different 
models of the l/f noise in Section 3.2.2. 
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3.2.1 PC A results 

The instrumental noise is simulated as explained in Section 2.4, and 
we apply the PCA method to the maps of the instrumental noise. 
The 1 // noise is computed with the knee frequency /knee = 1 mHz, 
which is thought to be achievable using balanced correlator re¬ 
ceivers (Jarosik et al. 2003; Bersanelli et al. 2010). Note that a dif¬ 
ferent scanning strategy than a drift-scan strategy can remove the 
l/f noise in the map-making, but for a transit telescope, such as 
BINGO, we will rely on component separation and the smoothness 
of the frequency spectrum. In the top panel of Fig. 4 we plot the 
power spectra of the noise maps with different number of modes 
removed: 1, 2 and 3. The spectra of the maps are computed to 
i <1000 using PolSpice (Szapudi, Prunet & Colombi 2001; Chon 
et al. 2004; Challinor & Chon 2005). This code computes corre¬ 
lation functions and estimates the power spectra by integrating the 
resampled correlation function using Legendre-Gauss integration. 
The power spectra are corrected for the effect of the cut sky and 
for the beam and pixel window functions. In order to remove the 
ringing in the power spectra, we apodize the correlation function 
with a Gaussian of width 15°. Subtracting one mode does not re¬ 
move the 1 //noise sufficiently well, but the thermal noise level can 
be reached by removing 2 modes as displayed in Fig. 4. This plot 
shows the residuals between the recovered thermal noise and the 
input thermal noise. For the case of two removed modes, the resid¬ 
ual is significantly lower than the input thermal noise at all scales. 
It shows that we can recover the thermal noise model sufficiently 
well using principal component analysis by subtracting at least 2 
principal modes. 

3.2.2 Results with different noise models 

In order to test the efficiency of the PCA method, we compute the 
instrumental noise with different values of the knee frequency /knee 
between 1 mHz and 10 Hz. Note that a knee frequency of 1 mHz 
might be expected for a pseudo-correlation receiver and that 10 Hz 
is a worst-case scenario for a single channel radiometer. We find 
that the residual noise after removing 2 principal modes is indepen¬ 
dent of the input knee frequency and, hence, the PCA method is ro¬ 
bust. We emphasise that we have assumed a flat frequency spectrum 
for the 1// noise (i.e. it affects all frequency channels equally). The 
efficiency of the noise cleaning method depends on this assumption 
and on assuming a perfect calibration. We will quantify the success 
of this method as the function of the smoothness of the l/f noise 
in Section 4. 

3.3 Inclusion of foreground emission 

In this Section, we show the results for the detection of the HI signal 
from the total intensity maps in the presence of instrumental noise 
and foreground emission, using the parametric fitting and principal 
component analysis methods. We compare both methods in Sec¬ 
tion 3.3.1 and we focus on the results of the principal component 
analysis method in Section 3.3.2. We look at the residuals of the 
reconstructed cosmological signal in Section 3.3.3. 

3.3.1 Parametric fitting results 

Here we present the results of the parametric fitting method. The 
extraction of the HI signal is done using only the frequency infor¬ 
mation. The fit is made for each pixel. In Fig. 5, we show the mea¬ 
surements as a function of frequency, for a random line-of-sight. 



Multipole t 

Figure 4. Power spectra of the simulated noise maps after applying princi¬ 
pal component analysis. We show the input 1// noise {red solid line) and the 
input thermal noise {green solid line), the reconstructed noise after apply¬ 
ing principal component analysis with 1 mode removed {black dashed line) 
and with 2 modes removed {orange dashed line). We display the power 
spectrum of the difference between the maps of the input thermal noise and 
of the reconstructed thermal noise after PCA with 2 modes removed {blue 
dotted line). It is below the input thermal noise at all scales. 


with the synchrotron model 3 quantified in Table 1 and a back¬ 
ground of unresolved point sources {S < 100 mJy) as explained in 
Section 2.1.2. The result is averaged over 20 realisations of the 
instrumental noise (thermal and l/f noise). The top panel repre¬ 
sents the simulated measurements and the reconstructed foreground 
emission with parametric fitting, highlighting the smooth compo¬ 
nent of the foreground. The bottom panel represents the recovered 
cosmological signal with parametric fitting and principal compo¬ 
nent analysis after removing 7 modes compared to the input one. It 
shows that the parametric fitting, while superficially in agreement 
with the input signal, does not provide an accurate fit to the signal 
of interest compared to the PCA method. 

3.3.2 Results of PCA applied to sky maps 

Fig. 5 shows that PCA induces a small offset in the reconstructed 
HI signal: some cosmological signal leaks into the reconstructed 
foreground and noise components. However, with this method, the 
HI signal is well recovered with a relative error of ~7%. 

Here we will quantify the impact of the foreground residuals 
on the reconstructed HI signal after applying the PCA method to 
the simulated sky maps. In order to check the effectiveness of the 
foreground removal method, we define the relative error by 

|7hi(v) - THi(v)|/an, (23) 

where 7hi(v) is the true HI signal at frequency v, Thi(v) the re¬ 
covered signal, and an the standard deviation of the thermal noise. 
Fig. 6 represents the relative error as a function of the number of 
subtracted modes for a simulation with HI signal, l/f and thermal 
noise and different models of foreground. We apply the PCA tech¬ 
nique to each of the foreground models discussed in Section 2.1.1. 

Since the foreground and the l/f noise spectra do not con¬ 
tain sharp features, we can expect that they are well described by a 
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Frequency (MHz) 

Figure 5. Top panel: brightness temperature of one particular line-of-sight 
as a function of frequency composed of sky emissions and instrumental 
noise (blue solid line), averaged over 20 realisations and the fit based on 
a parametric fitting {red dashed line). Bottom panel: the recovered HI sig¬ 
nal from parametric fitting {red dashed line) and from principal component 
analysis {black dashed line). We show the input HI {blue solid line) and the 
thermal noise {green solid line). 


small number of eigenvectors, so the eigenvalues are much larger 
for the first few principal components. This implies that a small 
number of components contains almost all of the foreground emis¬ 
sion and the 1// noise. Fig. 6 shows a fall-off of the amplitude 
of the eigenvalues with an increase in the value of the number of 
principal components. This steep drop means that the spectra are 
dominated by relatively few components, which are related to the 
foreground and smooth instrumental contamination. Furthermore, 
this figure shows that the foreground model has an impact on the 
extraction of the HI signal as the most complex foreground model 
(model 3) requires more modes to be removed in order to subtract 
the same level of foreground contamination. With the same noise 
model, the first foreground model requires the removal of at least 3 
principal modes, the second model, 4 modes and the third model, 7 
modes. 

We notice that, when a larger number of modes is removed 
beyond a certain threshold, the relative error associated with each 
foreground model begins to increase. This can be understood in 
terms of the component separation method inducing a leakage of 
the cosmological signal in the foreground eigenvectors. This be¬ 
comes even more important when we remove a larger number of 
principal modes from the initial maps. For the most realistic fore¬ 
ground model, model 3, the cosmological signal is well recovered 
with a relative error of ^1% after subtracting 7 modes and the per¬ 
centage error increases to ^9.5% with 11 modes removed. 

To evaluate the performance of the PC A method, we calculate 


the amount of leakage of the HI signal into the subtracted modes 
using Eq. 20. The relation between the maps of the HI signal at each 
frequency channel Sm and the orthogonal matrix P which contains 
the principal modes of the maps is given by the eigenfunctions (|) 

(^ = P^5hi. (24) 

Finally, we obtain the maps of the leakage of the HI signal us¬ 
ing the orthogonal matrix Pc, which contains only the eigenvectors 
removed from the initial maps 

Sm = Pc^- (25) 

Removing 7 principal modes with the PCA method enables one to 
recover the HI signal. However, this induces a loss of the cosmo¬ 
logical signal of ~5%. In order to determine the foreground com¬ 
ponent that has the most impact on the foreground separation, we 
generate maps of the sky emission with only synchrotron emission 
and HI signal, and maps with a background of unresolved point 
sources and HI signal. The maps with synchrotron emission require 
the removal of 3 modes in order to extract the signal of interest, and 
we obtain a leakage of the HI signal in the removed principal modes 
of ~2.4%. To recover the HI signal from maps with only a back¬ 
ground of unresolved point sources, the subtraction of 1 mode is 
required and leads to a leakage of ~1.6% of the cosmological sig¬ 
nal. Subtracting too many principal modes induces a significative 
loss of the cosmological signal, while when an insufficient number 
of modes is removed, the recovered HI signal will be still affected 
by the foreground contamination. 

Up to now we have assumed a perfect smoothness of the fore¬ 
ground spectra. This assumption could be broken in the presence 
of instrumental systematic effects such as not sufficient knowledge 
of the beams, imperfect polarisation purity and mis-calibration will 
affect the results of the component separation methods, adding ad¬ 
ditional uncertainties. To measure the BAO wiggles from the HI 
power spectrum, we require that the statistical error dominates the 
errors from foreground cleaning methods and calibration. For the 
present analysis, we neglect calibration systematics, postponing 
discussion of the required calibration accuracy needs to be done 
to a future paper. We are however aware that accurate calibration of 
bandpasses and beam polar diagrams is essential. Quantifying cali¬ 
bration requirements will be done with an end-to-end pipeline. We 
are also investigating a more complex foreground cleaning method, 
which uses combined spatial and spectral filtering techniques based 
on the expertise from the CMB (e.g. Leach et al. 2008; Remazeilles, 
Delabrouille & Cardoso 2011 and Olivari et al. in prep.). 


3.3.3 Maps of the reconstructed HI signal and power spectra 

We perform the foreground cleaning using the PCA method with 7 
modes subtracted and the parametric fitting. In Fig. 7, we display 
the maps of the recovered HI signal for the foreground model 3. The 
simulations contain instrumental noise (1// and thermal noise). To 
compare with the results of the cleaning methods, we plot in the 
first strip the true HI signal. The second, third and fourth strips 
show the reconstructed HI signal after 1, 3 and 7 principal modes 
are removed respectively. One can notice the similarity between 
the reconstructed HI signal shown in the fourth strip and the strip 
of the true HI signal. These two maps show that we can extract 
the cosmological signal from a highly contaminated map. The fifth 
strip represents the recovered HI signal after applying parametric 
fitting. Noting that the temperature scales are different, we observe 
significant differences between the input HI strip and that recovered 
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Figure 6. The percentage relative error of the thermal noise as a function 
of the number of modes removed. We plot the error for a foreground model 
with a background of unresolved point sources and Galactic synchrotron 
emission with a constant index p = —2.8 (red dashed line), with the in¬ 
troduction of the random spatial distribution of the index p across the sky 
(green dashed line), and with correlated spatial and spectral variations of 
the index p plus a background of unresolved point sources (black dashed 
line). These models correspond respectively to models 1, 2 and 3 explained 
in Table 1. The relative error increases significatively above 10 modes due 
to the leakage of the cosmological signal into the reconstructed foreground 
and correlated noise components. 
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with parametric fitting. It is clear that parametric fitting does not 
provide a sufficiently good fit to foreground and noise components. 

It is evident that parametric fitting does not remove the fore¬ 
ground completely, thus resulting in foreground leaking into the 
reconstructed cosmological signal. To highlight the comparison be¬ 
tween parametric fitting and PCA, we show in Fig. 8 the dispersion, 
pixel by pixel, between the recovered HI signal as a function of the 
input HI signal, obtained from parametric fitting and from PCA. 
This plot shows that parametric fitting is much less effective than 
the PCA method and induces a bias in the reconstructed HI signal. 

We can quantify the leakage of the thermal noise and of the 
cosmological signal into the reconstructed foreground and noise 
components by calculating their power spectra and comparing them 
to their input power spectra. In the bottom panel of Fig. 9, we dis¬ 
play the power spectra of the true and the recovered HI signal after 
the removal of 7 modes. This figure shows that both parametric fit¬ 
ting and PCA methods remove several orders of magnitude of fore¬ 
ground contamination and 1// noise, but PCA gives lower residuals 
than the parametric fitting method. We plot the cosmological sig¬ 
nal leakage into the foreground and 1// noise reconstruction in the 
bottom panel. The power spectra of the thermal noise leakage and 
the HI signal leakage are lower than the input HI signal at all scales, 
thus, with PCA, it is feasible to extract the HI signal from a highly 
contaminated foreground map. 


4 REQUIREMENTS ON FOREGROUND AND 

INSTRUMENTAL NOISE FREQUENCY SPECTRA 

For almost all component separation techniques developed for in¬ 
tensity mapping data analysis, the efficiency of these methods de¬ 
pends on the spectral smoothness of the foreground and that of 
the 1// noise and the bandpass calibration. We are confident that 
the foregrounds have sufficiently smooth spectra in frequency. This 
characteristic enables us to remove the components correlated in 
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Input HI signal (mK) 


Figure 8. Scatter plot of the recovered HI signal from parametric fitting 
(blue markers) and from principal component analysis (red markers) as a 
function of the input HI signal. The parametric fitting is more noisy than 
the principal component analysis and also appears to be biased. The black 
dashed line represents the perfect correlation between the recovered HI sig¬ 
nal and the true signal. 


frequency and, therefore, to recover the HI signal with the in¬ 
strumental white noise. However, some structure in the receiver 
band-passes is inevitable caused by amongst other things, stand¬ 
ing waves, frequency variations in the receiver gain temperature 
and spectrally dependent beam patterns. The effect of these can 
all be reduced by careful calibration but it is important to know 
how good this calibration has to be. In this section we define some 
requirements on the smoothness of the instrumental bandpass. In 
the following, for the foreground emission, we consider only the 
brightest component; i.e. the Galactic synchrotron emission. 

To quantify the effect of a non-smooth bandpass, we add to 
the measured frequency spectrum a sinusoidal wave defined by 

(|)(v)=Asin(|^j, (26) 

where A is the amplitude, v the frequency of observation and Av 
the wavelength. We explore a range of A between 1 and 150 mK 
and a range of Av between 1 to 300 MHz. 

We show the modified spectra for different values of A and Av 
in Fig. 10. To highlight the impact of the addition of the sinusoidal 
wave, we divide the resulting spectrum by the original one. We 
see curvature and/or oscillations in the resulting spectra. A higher 
value of Av leads to a curvature of the spectrum, similar to a stand¬ 
ing wave, whereas a smaller value induces a sinusoidal wave that 
behaves in a similar way to noise, when Av is smaller than the fre¬ 
quency channel width. 

In what follows, we simulate the maps generated by the instru¬ 
ment with the same model of the Galactic synchrotron emission, HI 
signal and instrumental noise (thermal and 1// noise) and we add 
Eq. 26 the sinusoidal wave to the frequency spectrum of the gen¬ 
erated data. In order to extract the signal of interest, we apply the 
principal component analysis to the maps. In Fig. 11, we plot the 
relative error of the recovered HI signal as a function of the am¬ 
plitude A and Av after applying PCA with 6 modes subtracted. The 
colour bar represents the amplitude of the relative error between the 
recovered HI signal and the true signal. The smoothness of the fre¬ 
quency spectrum, i.e. the value of Av, has a significative impact on 
the efficiency of the cleaning methods. The relative error increases 
with a small value of Av, which corresponds to a sinusoidal wave 
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Figure 7. Five versions of a strip in declination of the HI signal. The input cosmological signal is shown in the first strip and the reconstructed signal from 
principal component analysis in strips 2, 3 and 4 (after 1, 3 and 7 modes removed respectively) and from parametric fitting in strip 5. Notice the different 
colour bar scales. 


with a period shorter than the frequency channel width. In order 
not to be affected by the variation of the bandpass, the value of Av 
has to be lower than 100 MHz and the amplitude A has to be below 
45 mK. With the values A < 40 mK and Av < 100 MHz, we find a 
relative error <13% after 6 modes are subtracted with the PCA. 
In absolute terms, after subtracting 6 principal modes, we obtain 
residuals lower than 0.1 mK, which means that the HI signal can be 
detected. 


Finally, we perform simulations varying the number of fre¬ 
quency channels used to perform the PCA. We consider 20 fre¬ 
quency channels (15 MHz channel bandwidth) and 200 frequency 
channels (1.5 MHz channel width) and we test different values 
of Av. We choose the amplitude of the sinusoidal wave to A = 
120mK. Fig. 12 shows the relative error between the recovered HI 
signal and the input signal as a function of the smoothness of the 
frequency spectrum Av after removing 6 and 7 principal modes. 
We find that the PCA method does better with a larger number of 
channels. The relative error is < 7% for 6 removed modes when 
we have 200 frequency channels for all values of Av between 1 
to 400 MHz. The reason for the improvement when more channels 
are added can be understood by the fact that the frequency band 
is better sampled. Thus, as long as we have a frequency spectrum 
with slow oscillations, or enough frequency channels to sample the 
spectrum with sufficient accuracy, the smoothness of the bandpass 
does not constitute an issue for the foreground and the 1// noise 
cleaning methods. An amplitude around 40 mK requires the band¬ 
pass to be calibrated to an accuracy of better than 1 part in 1000. 
However, one would expect to calibrate at least every day so we 
will only require a dynamic range of 1 part in 50. 


5 CONCLUSION 

In Section 2, we have simulated the kind of data which would 
be produced in a single-dish intensity mapping experiment like 
BINGO. We have adopted a sky model with Galactic synchrotron 
emission and a background of point sources (< lOmJy). We have 
added both white and 1// noise. This simulation can be generalised 
to any intensity mapping experiment with a drift scan strategy. This 
work has yielded several results concerning the challenge of any 
future single-dish experiment in terms of foreground, atmospheric 
and instrumental noise. 

We have made some estimates of the total amplitude of the 
atmospheric noise and of its fiuctuations arising from turbulences 
in emission of the water vapour. At 1 GHz, we have found that the 
r.m.s. noise level induced by these fiuctuations is ~ 0.01 mK below 
the instrumental noise level for the BINGO experiment. Thus, for 
an observing frequency around 1 GHz, the atmospheric contamina¬ 
tion should not constitute a challenge for a single-dish experiment. 

We have investigated the problem of cleaning foreground con¬ 
tamination in order to recover the HI signal. We have focused on 
simple methods based on the spectral information contained in 
the signal. We have considered two different foreground cleaning 
methods, a non-blind method with parametric fitting and a blind 
method with PCA. This last method does not require any assump¬ 
tion on the physics of the foreground emissions. We have shown 
that the parametric fitting method does not provide an accurate fit 
to the data for a simulation with instrumental noise and foreground 
emissions. In contrast, by subtracting 7 principal modes, with a re¬ 
alistic foreground model, PCA enables to reach the HI signal level. 
Thus, on simulated data, the application of the PCA method shows 
that it is feasible to extract the cosmological signal across a wide 
range of multipoles and redshifts. 

However, we have found that PCA induces a small offset in 
the reconstructed the HI signal of ^5%. This result is confirmed by 
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Figure 9. Top panel: the power spectra of the simulated maps after apply¬ 
ing foreground cleaning. This plot shows the simulated HI signal (blue solid 
line), the thermal noise (green solid line) and the sky emission (cyan solid 
line). We plot the results of the foreground cleaning from principal com¬ 
ponent analysis (black dashed line) and from parametric fitting (red dashed 
line). Bottom panel: the power spectra of the leakage of the HI signal (dotted 
light blue line) and of the noise (gold dashed line) after applying PCA. The 
principal component analysis method clearly makes it possible to extract 
the HI signal. 



Figure 10. The corrupted spectrum divided by the original, undistorted, 
spectrum for different versions of the sinusoidal wave. A small value of 
Av induces a spectrum that fluctuates in frequency in a similar way to ran¬ 
dom noise, while a large value of Av leads to a curvature of the frequency 
spectrum. 
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Figure 11. Relative error as a function of the amplitude A and the period Av 
of the sinusoidal wave after foreground and noise subtraction after applying 
principal component analysis (6 modes removed). The values A and Av of 
the sinusoidal waves are indicated on the axes of the plot. The colour bar 
gives the percentage error relative to the noise. 
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Figure 12. The percentage relative error as a function of the period Av of 
the sinusoidal wave after foreground removal. We show the relative errors 
for 20 frequency channels with 6 modes removed (blue markers), for 20 
frequency channels with 7 modes removed (green markers) and for 200 
frequency channels with 6 and 7 modes removed (red and black markers) 
respectively. 


other recent works (Alonso et al. 2015), which have shown that it 
is possible to recover successfully HI signal from a highly contam¬ 
inated foreground map with a blind method but that the foreground 
removal induces a bias in the HI power spectrum. In spite of this 
HI leakage into the removed contamination, a modification of the 
global shape of the power spectrum should not modify the posi¬ 
tions of the BAOs wiggles (Wolz et al. 2014) and therefore will not 
have significant impact on our ability to measure the BAO scale. 
If the PCA method is chosen, we can correct for this small bias 
by calibrating it with simulations. We are also investigating a more 
complex cleaning method, which utilises both the spectral and spa¬ 
tial information (Olivari et al., in prep.), and therefore should be 
more accurate. 

Given the assumption that the frequency spectrum of the 1// 
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noise is flat, or at least only has slow variations across the band, 
we have shown that PCA can cope with 1// noise with relatively 
high knee frequencies. It is easily removed from the noise maps 
by subtracting 2 principal modes and the 1// noise contamination 
is dampened down to the level of the thermal noise. This result 
implies that it might be possible to simplify the design of an exper¬ 
iment, such as BINGO, by removing the correlation receivers. In 
subsequent work, we plan to investigate more complex 1// noise 
models. 

The effectiveness of the blind cleaning methods depends on 
the spectral smoothness of the foregrounds and the instrumental 
1// noise. In order to challenge our cleaning methods, we have 
considered instrumentally distorted frequency spectra of the com¬ 
ponents we want to remove. We have shown that a sinusoidal 
distortion of the frequency spectrum can induce higher residuals 
on the maps. However, for a sinusoidal wave with an amplitude 
A < 45 mK and period Av < 100 MHz, after subtracting six princi¬ 
pal modes from the initial maps, the residuals of the noise and the 
sky emission are below 0.1 mK, which could enable the successful 
extraction of the HI signal. Still, accurate calibration will be criti¬ 
cal for the success of measuring B ary on Acoustic Oscillations and 
probing the expansion history of the Universe. How good this cali¬ 
bration needs to be will be investigated in a future paper. The results 
of this paper show that the principal component analysis method is 
a promising tool for the extraction of the cosmological signal in any 
future HI intensity mapping experiment. 
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APPENDIX A: ATMOSPHERIC NOISE 

We quantify the total contribution of the atmospheric noise in Sec¬ 
tion A1 and we discuss the contribution of the atmospheric turbu¬ 
lences induced by the water vapor inhomogeneities in Section A2. 


The total contribution of the atmospheric emission Tem to the 
power received on one receiver is related to the optical depth by 

7;m = Tatm(l-^“"). (A2) 

The atmosphere leads to an attenuation of the signal from a source 
with a fiux density 5'em- For an average optical depth x, the fiux 
received by the instrument is 

‘^abs ~ • (A3) 

Eor T << 1, we obtain S'abs 5'em(l — "c). Thus, the attenuation of 
the source by the atmosphere is at a level of one percent. Eor this 
reason, we can consider only the atmospheric emission, which is 
much larger than the absorption. 

The atmospheric emission contributes to the system tempera¬ 
ture of the instrument, adding to the noise from the emissions of 
the instrument subsystems and from the ground pick-up and sky 
(astronomical signal). In the case of a small optical depth t ^ 0.05, 
from Eq. A2, we can write Tem ~ Tatm'C. The theoretical emission 
is given by 

7’em(lGHz) = 1.8+ 3.9 X 10"^PWV (mm). (A4) 

In days with favourable atmospheric conditions PWV < 2 mm, 
we find Tern = 1-81 K, and in days with bad weather PWV > 
4 mm, Tern — 1 .82 K. Prom the models of the atmospheric emission 
(Danese & Partridge 1989), Staggs et al. (1996) predicts a similar 
contribution with Tem = 1.7 K. It varies so little because T^v '^ox 
when the weather is stable. 

A2 Atmospheric fluctuations 

The fluctuating atmospheric signal in the timelines depends on the 
precipitable water vapour level, the turbulence height and the thick¬ 
ness and the wind speed of the layer. According to the Kolmogorov 
model of turbulence (Tatarskii 1961), we can assume that the fluc¬ 
tuations in water vapour occur in a turbulent layer at a height hav 
with a thickness 5/z, with hav » ^h. The resulting power spectrum 
from the turbulence in the timelines is assumed to be a power-law 
with index b = —11/3 on small scales and Z? = — 8/3 at large scales, 
and is given for both cases by (Lay & Halverson 2000) 


Al Total atmospheric absorption and emission 

Eor our analysis of the effects of absorption and emission below, 
we consider a layer of the water vapour at a temperature Tatm and 
an optical depth x at a given observing frequency. Around 1 GHz, 
the optical depth is dominated by oxygen with a small contribution 
from water, usually quantified with the precipitable water vapour 
(PWV). PWV is defined as the height of liquid water if all the water 
vapor would be condensed in a column of the atmosphere. Hence, 
it can undergo changes on a daily basis and between different sites. 
The total zenith optical depth can be calculated by the following 
relation at one frequency v (Butler 2002) 

X = 10-3 X (Xox + XwvPWV (mm)), (Al) 

where Xqx is the oxygen optical depth and Xwv the water vapor op¬ 
tical depth. These parameters are constant and depend on the fre¬ 
quency. We use the ATM software (Pardo, Cernicharo & Serabyn 
2001) to determine their values for an altitude of the site at 200 m, 
a ground level pressure of ^1000 mbar, a water vapour scale height 
of 2 km and a ground level temperature Tatm = 260 K, and we find 
Xqx = 7.1, Xwv = 0.015. This is consistent with other predictions in 
the literature (Smith 1982). 


< T^{ax,ay) >= 


A 

sine 



if- 


/la 


f 2 

<(^aj+a;j <ainner, 


2Ah sine 


(A5) 


< T^{ax,ay) >= 



/ 7 9\*/2 

if OCouter ^ “1“ ^>7 J ^ 


/lav 

2Ah sine’ 


(A6) 


where ainner and aouter are the inner and outer scales of the turbu¬ 
lences, Ah is the outer scale size, and e is the elevation angle. The 
coefficients A and A are the amplitudes of the turbulence. The pa¬ 
rameters axj are the angular wavenumbers, related to the spatial 
wavenumbers kxj by ax^y = ^jc^/iav/sine. 

In order to predict the residual level of the brightness tempera¬ 
ture fluctuations for a single dish at sea level, we follow the Church 
model (Church 1995). We assume that the effective beam area is 
described by a Gaussian beam 


71W^ (z) 


2(x^±yf] 

w2(z) J ’ 


(A7) 
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where w(z) = wq (^1 + ■ 

The parameter wq is the beam waist radius related to the resolution 
0FWHM by wo = Xy/2\og2/ (ttGfwhm)- 

In the limit of a long averaging time, the Church model 
(Church 1995) gives the rms fluctuations in antenna temperature 
for a single-dish according to the outer scale Lq by 



where Tatm is the atmosphere temperature, L = 3Lo, Zu is the height 
of the atmosphere and z\ the breakdown of the assumption z 
v,y. The parameter Ca is a measurement of the amplitude of the 
turbulent fluctuations and is given at 15 GHz by 

C^ = 2.0xl0"''^m"*/^exp^-^^. (A9) 

To simplify the equations, we do not take into account the effect 
of the wind speed as done in Lay & Halverson (2000). This adds a 
time dependence, which allows the measured fluctuations to de¬ 
crease with larger integration times. The brightness temperature 
of the fluctuations induced by the water vapour molecules can be 
scaled to other frequencies, based on the emissivity spectrum of 
water vapour. The water vapour spectrum is composed of spec¬ 
tral emission lines overlaid on a background emission, which in¬ 
creases as the square of the frequency (Smith 1982), so we scale 
Eq. A8 as our frequency of observation to the fourth power. We as¬ 
sume a model of the telescope with Gfwhm = 40 arcmin at 1 GHz, 
a temperature of the atmosphere of Tatm = 270 K and a maximum 
height of 2 km. We use the value suggested in Church (1995) for 
Lq = 10 m. With these assumptions, we And ATatm ^ 0 01 rnK. 
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